Fisheries induced shifts in sea snake assemblages along the Konkan coast of India
Supplementary materials B: Analysis and Figures
options(max.print="75")
knitr::opts_chunk$set(echo = T, error = T, warning = F, message = F, tidy = T, cache = T,
dpi = 300, fig.align = "center", eval = T)
#Required libraries
suppressPackageStartupMessages(library(plyr))
suppressPackageStartupMessages(library(tidyverse))
suppressPackageStartupMessages(library(lubridate))
suppressPackageStartupMessages(library(RColorBrewer))
#Importing sea snake bycatch data
snakes = read.csv("./Data/Rao et al_data_sea snake bycatch.csv")
# Fixing tow duration for 2018 - 2019
snakes <- snakes%>%
mutate(Month = factor(months(dmy(Date), abbreviate = T), levels = month.abb))%>%
group_by(Year, Gear.Type)%>%
mutate(n.hauls = ifelse(is.na(No..of.Hauls), #filling in missing data
mean(No..of.Hauls, na.rm = T), No..of.Hauls),
haul.time = ifelse(is.na(Average.Haul.Duration..Hours.),
mean(Average.Haul.Duration..Hours., na.rm = T),
Average.Haul.Duration..Hours.))%>%
ungroup()%>%
mutate(Tow.duration.hours =
ifelse(is.na(Tow.duration.hours), # only for TD not recorded
ifelse(Gear.Type == "GillNet", #for gillnets
(haul.time+1.75), #soak time correction
ifelse(Gear.Type == "Trawler", #for trawlers
n.hauls*haul.time, # total tow duration
Tow.duration.hours) #return unchanged
), #return unchanged
Tow.duration.hours)# return same for TD recorded
)
#setting ggplot theme
theme_set(theme_bw()+
theme(legend.position = "top",
text = element_text(size = 15)
)
)
# function calculating mcfadden's pseudo r2 for glms
mcfadden <- function(x){
r2 <- 1 - (x$deviance/x$null.deviance)
}
# function for calculating effect size for Z - test
proptest.r <- function(x, y){
r <- x$statistic/sum(y$N)
}Composition of sea snakes in bycatch
Sampling effort
We calculated the number of trips sampled and fishing effort per trip (Tow duration in hours) to determine sampling effort per month of each year of sampling.
#Number of boats sampled accross years
snakes%>%
group_by(Year, Month)%>%
distinct(Date, Boat.Name, .keep_all = T)%>%
#filter(!Boat.Name == "")%>%
count(name = "n.trips")%>%
spread(Month, n.trips)| Year | Jan | Feb | Mar | Apr | May | Oct | Nov | Dec |
|---|---|---|---|---|---|---|---|---|
| 2016 | 2 | 8 | 8 | 3 | 6 | NA | 4 | 10 |
| 2017 | 12 | 19 | 15 | NA | NA | NA | 1 | 6 |
| 2018 | 91 | 56 | 59 | 15 | 22 | 18 | 41 | 44 |
# Fishing effort accross years
snakes%>%
filter(Gear.Type == "GillNet" | Gear.Type == "Trawler")%>%
group_by(Gear.Type, Year, Month)%>%
distinct(Date, Boat.Name, .keep_all = T)%>%
summarise(sampling.effort = round(sum(Tow.duration.hours, na.rm = T), 2))%>%
spread(Month, sampling.effort, fill = c(0))| Gear.Type | Year | Jan | Feb | Mar | Apr | May | Oct | Nov | Dec |
|---|---|---|---|---|---|---|---|---|---|
| GillNet | 2016 | 9.00 | 29.00 | 24.00 | 7.00 | 21.00 | 0.00 | 10.00 | 36.50 |
| GillNet | 2017 | 3.00 | 27.00 | 9.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 |
| GillNet | 2018 | 206.91 | 112.87 | 128.63 | 25.88 | 51.63 | 43.51 | 62.26 | 109.38 |
| Trawler | 2016 | 0.00 | 8.00 | 9.00 | 9.00 | 8.00 | 0.00 | 0.00 | 0.00 |
| Trawler | 2017 | 84.00 | 105.00 | 52.50 | 0.00 | 0.00 | 0.00 | 0.00 | 15.00 |
| Trawler | 2018 | 49.00 | 50.91 | 30.50 | 11.75 | 0.00 | 0.00 | 11.00 | 5.50 |
snakes%>%
filter(Gear.Type == "Trawler" | Gear.Type == "GillNet")%>%
group_by(Gear.Type)%>%
distinct(Date, Boat.Name, .keep_all = T)%>%
summarise(n.trips = n(),
sampling.effort = round(sum(Tow.duration.hours, na.rm = T), 2),
n.boats = length(unique(Boat.Name)))| Gear.Type | n.trips | sampling.effort | n.boats |
|---|---|---|---|
| GillNet | 329 | 916.57 | 136 |
| Trawler | 65 | 449.16 | 6 |
Sampling was not conducted trip wise for 2016 and most of 2017. Hence, effort data is unavailable.
Assemblage composition
We computed the abundance of each species of snake encountered over the whole sampling duration and across years
snakes%>%
group_by(Species, Year)%>%
count(name = "Abundance")%>%
spread(Year, Abundance)| Species | 2016 | 2017 | 2018 |
|---|---|---|---|
| Acrochordus granulatus | NA | 1 | NA |
| Hydrophis curtus | 75 | 38 | 123 |
| Hydrophis cyanocinctus | NA | NA | 3 |
| Hydrophis schistosus | 178 | 215 | 521 |
| Hydrophis viperinus | NA | 1 | 1 |
Varition in catch per unity effort (CPUE) of sea snakes across years in each gear
We considered only the two dominant species for in depth analysis. Catch Per Unit Effort (CPUE) for each species in each sampled trip was calculated as ‘Abundance of snakes encountered over Total Tow Duration or Soak Time’. We tested the effect of Tow duration on CPUE of each species using a log – linear regression. We compared the catch rates (CPUE) of each species across and within gears using log – linear regression.
CPUE.grsp <- snakes%>%
filter(Gear.Type == "GillNet" | Gear.Type == "Trawler")%>%
filter(Species == "Hydrophis curtus" |
Species == "Hydrophis schistosus",
Boat.Name != "")%>% # keeping dominant species
droplevels()%>%
group_by(Year, Gear.Type, Date, Boat.Name, Species)%>%
summarise(Abn = n(), # abundance of sea snakes per trip
fe = mean(Tow.duration.hours, na.rm = T))%>% # fishing effort, mean for 2016 - 17
complete(nesting(Year, Gear.Type, Date, Boat.Name, fe),
Species, fill = list(Abn = 0))%>% # completing species x trip combinations
mutate(cpue = Abn/fe)%>%
group_by(Boat.Name)%>%
mutate(n = n())%>%
ungroup()
CPUE.grsp%>%
group_by(Gear.Type, Species)%>%
summarise(Mean.CPUE = mean(cpue, na.rm = T))%>%
spread(Species, Mean.CPUE)| Gear.Type | Hydrophis curtus | Hydrophis schistosus |
|---|---|---|
| GillNet | 0.4243868 | 0.4647042 |
| Trawler | 1.1981818 | 1.7018479 |
CPUE.grsp%>%
group_by(Gear.Type,Species)%>%
nest()%>%
mutate(shapiro = map(data, ~shapiro.test(.$cpue)),
sumry = map(shapiro, broom::tidy))%>%
select(sumry)%>%
unnest()| Gear.Type | Species | statistic | p.value | method |
|---|---|---|---|---|
| GillNet | Hydrophis schistosus | 0.4785970 | 0.0000000 | Shapiro-Wilk normality test |
| GillNet | Hydrophis curtus | 0.7492041 | 0.0000001 | Shapiro-Wilk normality test |
| Trawler | Hydrophis curtus | 0.6313492 | 0.0001311 | Shapiro-Wilk normality test |
| Trawler | Hydrophis schistosus | 0.7890101 | 0.0003384 | Shapiro-Wilk normality test |
# testing difference in CPUE between trawlers and gillnets
library(lme4)
library(broom.mixed)
CPUE.grsp%>%
group_by(Species)%>%
nest()%>%
mutate(mod = map(data, ~lmer(cpue ~ Gear.Type + (1|Boat.Name), data = .)),
sumry = map(mod, broom.mixed::tidy))%>%
select(sumry)%>%
unnest()| Species | effect | group | term | estimate | std.error | statistic |
|---|---|---|---|---|---|---|
| Hydrophis schistosus | fixed | NA | (Intercept) | 0.4640419 | 0.0435137 | 10.664273 |
| Hydrophis schistosus | fixed | NA | Gear.TypeTrawler | 1.2257846 | 0.1431319 | 8.564023 |
| Hydrophis schistosus | ran_pars | Boat.Name | sd__(Intercept) | 0.0555275 | NA | NA |
| Hydrophis schistosus | ran_pars | Residual | sd__Observation | 0.6191548 | NA | NA |
| Hydrophis curtus | fixed | NA | (Intercept) | 0.4206151 | 0.1105983 | 3.803089 |
| Hydrophis curtus | fixed | NA | Gear.TypeTrawler | 1.0609298 | 0.3247348 | 3.267065 |
| Hydrophis curtus | ran_pars | Boat.Name | sd__(Intercept) | 0.4440595 | NA | NA |
| Hydrophis curtus | ran_pars | Residual | sd__Observation | 0.5613696 | NA | NA |
CPUE.grsp%>%
group_by(Species)%>%
nest()%>%
mutate(mod = map(data, ~lmer(cpue ~ Gear.Type + (1|Boat.Name), data = .)),
anva = map(mod, car::Anova),
sumry = map(anva, broom::tidy))%>%
select(sumry)%>%
unnest()| Species | term | statistic | df | p.value |
|---|---|---|---|---|
| Hydrophis schistosus | Gear.Type | 73.34249 | 1 | 0.0000000 |
| Hydrophis curtus | Gear.Type | 10.67372 | 1 | 0.0010867 |
# testing difference in CPUE between Species
CPUE.grsp%>%
group_by(Gear.Type)%>%
nest()%>%
mutate(mod = map(data, ~lmer(cpue ~ Species + (1|Boat.Name), data = .)),
sumry = map(mod, broom::tidy))%>%
select(sumry)%>%
unnest()| Gear.Type | effect | group | term | estimate | std.error | statistic |
|---|---|---|---|---|---|---|
| GillNet | fixed | NA | (Intercept) | 0.4244221 | 0.0437981 | 9.6904140 |
| GillNet | fixed | NA | SpeciesHydrophis schistosus | 0.0402291 | 0.0485110 | 0.8292770 |
| GillNet | ran_pars | Boat.Name | sd__(Intercept) | 0.0082453 | NA | NA |
| GillNet | ran_pars | Residual | sd__Observation | 0.3000941 | NA | NA |
| Trawler | fixed | NA | (Intercept) | 1.6265715 | 0.7494314 | 2.1704076 |
| Trawler | fixed | NA | SpeciesHydrophis schistosus | -0.3410434 | 0.6845323 | -0.4982138 |
| Trawler | ran_pars | Boat.Name | sd__(Intercept) | 1.1674491 | NA | NA |
| Trawler | ran_pars | Residual | sd__Observation | 1.5111595 | NA | NA |
CPUE.grsp%>%
group_by(Gear.Type)%>%
nest()%>%
mutate(mod = map(data, ~lmer(cpue ~ Species + (1|Boat.Name), data = .)),
anva = map(mod, car::Anova),
sumry = map(anva, broom::tidy))%>%
select(sumry)%>%
unnest()| Gear.Type | term | statistic | df | p.value |
|---|---|---|---|---|
| GillNet | Species | 0.6877004 | 1 | 0.4069477 |
| Trawler | Species | 0.2482170 | 1 | 0.6183334 |
#Plotting
snakes%>%
filter(Gear.Type == "GillNet" | Gear.Type == "Trawler")%>%
filter(Species == "Hydrophis curtus" |
Species == "Hydrophis schistosus")%>% # keeping dominant species
droplevels()%>%
group_by(Year, Gear.Type, Date, Boat.Name, Species)%>%
summarise(Abn = n(), # abundance of sea snakes per trip
fe = mean(Tow.duration.hours, na.rm = T))%>% # fishing effort, mean for 2016 - 17
complete(nesting(Year, Gear.Type, Date, Boat.Name, fe),
Species, fill = list(Abn = 0))%>% # completing species x trip combinations
group_by(Year, Gear.Type, Species)%>%
summarise(CPUE = mean(Abn/fe, na.rm = T), # Mean CPUE per trip
sd = sd(Abn/fe, na.rm = T))%>%
ggplot(aes(Year, CPUE, shape = Species, group = Species))+
geom_point(size = 3)+
geom_line(aes(linetype = Species), size = 1)+
scale_x_continuous(breaks = seq(2016, 2018, 1))+
labs(y = "Mean CPUE")+
facet_wrap(~Gear.Type)+
theme(legend.text = element_text(face = "italic"))
Seasonal variation in CPUE
We tested for seasonality in bycatch trends for each species in each gear using log – linear regressions with months as predictor variables.
# testing seasonality in CPUE by gear
CPUE.month <- snakes%>%
filter(Gear.Type == "GillNet" | Gear.Type == "Trawler")%>%
filter(Species == "Hydrophis curtus" |
Species == "Hydrophis schistosus")%>% # keeping dominant species
droplevels()%>%
group_by(Year, Month, Gear.Type, Date, Boat.Name, Species)%>%
summarise(Abn = n(), # abundance of sea snakes per trip
fe = mean(Tow.duration.hours, na.rm = T))%>% # fishing effort, mean for 2016 - 17
complete(nesting(Year, Month, Gear.Type, Date, Boat.Name, fe),
Species, fill = list(Abn = 0))%>% # completing species x trip combinations
mutate(CPUE = Abn/fe) # Mean CPUE per trip
CPUE.month%>%
ungroup()%>%
group_by(Gear.Type, Species)%>%
nest()%>%
mutate(mod = map(data, ~lmer(CPUE ~ Month + (1|Boat.Name), data = .)),
sumry = map(mod, broom::tidy),
tst = map(mod, car::Anova))%>%
select(sumry)%>%
unnest()%>%
select(term, estimate)%>%
spread(term, estimate, drop = T)| Gear.Type | Species | (Intercept) | MonthApr | MonthDec | MonthFeb | MonthMar | MonthMay | MonthNov | MonthOct | sd__(Intercept) | sd__Observation |
|---|---|---|---|---|---|---|---|---|---|---|---|
| GillNet | Hydrophis curtus | 0.3474585 | 1.2261151 | 0.1051078 | 0.0721653 | 1.6615975 | NA | 0.0344719 | 0.0343597 | 0.0000000 | 0.8967478 |
| GillNet | Hydrophis schistosus | 0.4313831 | 0.0670926 | 0.0095525 | 0.5993502 | 0.0363549 | 0.0028650 | -0.0718391 | -0.0360249 | 0.0949908 | 0.5935458 |
| Trawler | Hydrophis curtus | 0.3366162 | -0.1013220 | 0.4050505 | -0.0995791 | -0.1143939 | NA | 3.8300505 | NA | 0.0000000 | 0.6425533 |
| Trawler | Hydrophis schistosus | 0.7202427 | 1.6230328 | 0.9464240 | 0.3707643 | 0.4578930 | 0.1702764 | 2.0797573 | NA | 0.4216842 | 1.1841602 |
CPUE.month%>%
ungroup()%>%
group_by(Gear.Type, Species)%>%
nest()%>%
mutate(mod = map(data, ~lmer(CPUE ~ Month + (1|Boat.Name), data = .)),
sumry = map(mod, broom::tidy),
tst = map(mod, car::Anova))%>%
select(tst)%>%
unnest()| Gear.Type | Species | Chisq | Df | Pr(>Chisq) |
|---|---|---|---|---|
| GillNet | Hydrophis schistosus | 30.321943 | 7 | 0.0000829 |
| Trawler | Hydrophis schistosus | 8.054082 | 6 | 0.2341679 |
| GillNet | Hydrophis curtus | 23.811045 | 6 | 0.0005657 |
| Trawler | Hydrophis curtus | 67.330098 | 5 | 0.0000000 |
#plotting
snakes%>%
filter(Gear.Type == "GillNet" | Gear.Type == "Trawler")%>%
filter(Species == "Hydrophis curtus" |
Species == "Hydrophis schistosus")%>% # keeping dominant species
droplevels()%>%
mutate(Month = factor(Month, levels = month.abb))%>%
group_by(Year, Month, Gear.Type, Date, Boat.Name, Species)%>%
summarise(Abn = n(), # abundance of sea snakes per trip
fe = mean(Tow.duration.hours, na.rm = T))%>% # fishing effort, mean for 2016 - 17
complete(nesting(Year, Gear.Type, Date, Boat.Name, fe), Month,
Species, fill = list(Abn = 0))%>% # completing species x trip combinations
group_by(Month, Gear.Type, Species)%>%
summarise(CPUE = mean(Abn/fe, na.rm = T), # Mean CPUE per trip
sd = sd(Abn/fe, na.rm = T))%>%
ggplot(aes(Month, CPUE, shape = Species, group = Species))+
geom_point(size = 3)+
geom_line(aes(linetype = Species), size = 1)+
scale_y_log10()+
geom_vline(xintercept = "Jun", size = 1, linetype = "dotted")+
geom_vline(xintercept = "Aug", size = 1, linetype = "dotted")+
geom_text(aes(x = "Jul", y = 0.3, label = "Monsoon Ban"))+
labs(y = "Mean CPUE (log)")+
facet_wrap(~Gear.Type, ncol = 1)+
theme(legend.text = element_text(face = "italic"))
Length distribution in bycatch
We determined selectivity of gears with respect to SVL of individuals using a t – test for each species in each gear.
#plotting
snakes%>%
filter(Gear.Type == "GillNet" | Gear.Type == "Trawler")%>%
filter(Species == "Hydrophis curtus" |
Species == "Hydrophis schistosus")%>% # keeping dominant species
ggplot(aes(Snout.to.Vent..cm., fill = Gear.Type))+
geom_density(aes(linetype = Gear.Type),alpha = 0.3, size = 1)+
labs(x = "Snout to vent length (cm)", y = "Kernel Density")+
facet_wrap(~Species)+
theme(strip.text = element_text(face = "italic"))
# summary SVL
snakes%>%
filter(Species == "Hydrophis curtus" | Species == "Hydrophis schistosus")%>%
group_by(Species)%>%
skimr::skim(Snout.to.Vent..cm.)%>%
skimr::yank("numeric")Variable type: numeric
| skim_variable | Species | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
|---|---|---|---|---|---|---|---|---|---|---|---|
| Snout.to.Vent..cm. | Hydrophis curtus | 89 | 0.62 | 58.07 | 12.77 | 34.5 | 49.0 | 55.5 | 63.75 | 105.0 | ▃▇▃▂▁ |
| Snout.to.Vent..cm. | Hydrophis schistosus | 203 | 0.78 | 74.36 | 15.01 | 1.0 | 64.5 | 73.5 | 84.85 | 134.2 | ▁▁▇▃▁ |
#Z - test for proportion
snakes%>%
filter(Gear.Type == "GillNet" | Gear.Type == "Trawler")%>%
filter(Species == "Hydrophis curtus" |
Species == "Hydrophis schistosus")%>%
filter(!is.na(Snout.to.Vent..cm.))%>%
droplevels()%>%
select(Gear.Type, Species, Snout.to.Vent..cm.)%>%
group_by(Species)%>%
nest()%>%
mutate(t.test = map(data, ~t.test(Snout.to.Vent..cm. ~ Gear.Type, data = .)),
sumry = map(t.test, broom::tidy),
d = map(data, ~lsr::cohensD(Snout.to.Vent..cm. ~ Gear.Type, data = .)))%>%
unnest(sumry, d)%>%
select(estimate:p.value, d)| Species | estimate | estimate1 | estimate2 | statistic | p.value | d |
|---|---|---|---|---|---|---|
| Hydrophis curtus | -0.6316698 | 57.60896 | 58.24063 | -0.2802908 | 0.7797036 | 0.0489142 |
| Hydrophis schistosus | -1.6728765 | 74.32615 | 75.99903 | -1.4689341 | 0.1423528 | 0.1155450 |
Sex distibution in bycatch
We determined selectivity in terms of sex individuals using a two sample Z test for the proportion of females caught.
#Plotting
snakes%>%
filter(Sex == "Male" | Sex == "Female")%>% # sex recorded
filter(Gear.Type == "GillNet" | Gear.Type == "Trawler")%>%
filter(Species == "Hydrophis curtus" |
Species == "Hydrophis schistosus")%>% # keeping dominant species
ggplot(aes(Species, fill = Sex))+
geom_bar(position = "dodge", col = "black", width = 0.25)+
scale_fill_brewer(palette = "Greys")+
facet_wrap(~Gear.Type)+
theme(axis.text.x = element_text(face = "italic"))
# Z - test for proportion
snakes%>%
filter(Sex == "Male" | Sex == "Female")%>% # sex recorded
filter(Gear.Type == "GillNet" | Gear.Type == "Trawler")%>%
filter(Species == "Hydrophis curtus" |
Species == "Hydrophis schistosus")%>%
group_by(Gear.Type, Species, Sex)%>%
count()%>%
spread(Sex, n)%>%
mutate(N = Male+Female)%>%
group_by(Species)%>%
nest()%>%
mutate(prop = map(data, ~prop.test(.$Female, .$N)),
sumry = map(prop, broom::tidy),
r = map2(prop, data, ~proptest.r(x = .x, y = .y)))%>%
unnest(sumry, r)%>%
select(estimate1, estimate2, p.value, conf.low, conf.high, r)| Species | estimate1 | estimate2 | p.value | conf.low | conf.high | r |
|---|---|---|---|---|---|---|
| Hydrophis curtus | 0.5940594 | 0.2962963 | 0.0007665 | 0.1286188 | 0.4669074 | 0.0730372 |
| Hydrophis schistosus | 0.5082508 | 0.4444444 | 0.1848910 | -0.0282980 | 0.1559108 | 0.0034468 |
Significantly more male H. curtus caught in trawler than gillnets
Effect of effort on composition of sea snakes
We tested the effect of Tow duration on CPUE of each species using a log – linear regression.
#Extracting CPUE data from bycatch database
CPUE = snakes%>%
filter(Gear.Type == "GillNet" | Gear.Type == "Trawler")%>%
filter(!Boat.Name == "")%>% # removing points with no effort
filter(Species == "Hydrophis curtus" |
Species == "Hydrophis schistosus")%>% # keeping dominant species
droplevels()%>%
group_by(Gear.Type, Date, Boat.Name, Species)%>%
summarise(Abn = n(), # abundance of sea snakes per trip
fe = mean(Tow.duration.hours, na.rm = T))%>% # fishing effort, mean for 2016 - 17
complete(nesting(Gear.Type, Date, Boat.Name, fe),
Species, fill = list(Abn = 0))%>%
mutate(CPUE = Abn/fe,
log.Abn = log(Abn),
log.CPUE = log(CPUE),# log CPUE because CPUE is 0 inflated
log.fe = log(fe))# fe is 0 inflated
#Sample size
CPUE%>%
group_by(Gear.Type)%>%
distinct(Date, Boat.Name)%>%
count()| Gear.Type | n |
|---|---|
| GillNet | 255 |
| Trawler | 27 |
#Plot
CPUE%>%
ggplot(aes(log.fe, log.CPUE, shape = Species))+
geom_jitter(size = 3)+
geom_smooth(aes(linetype = Species), method = "lm", col = "black")+
scale_color_brewer(palette = "Greys")+
labs(y = "CPUE (log)", x = "Fishing Effort (log haul hours)")+
facet_wrap(~ Gear.Type, scales = "free_x")+
theme(legend.text = element_text(face = "italic"))
# testing
CPUE%>%
group_by(Gear.Type, Species)%>%
nest()%>%
mutate(mod = map(data, ~lmer(CPUE ~ fe + (1|Boat.Name), data = .)),
sumry = map(mod, broom::tidy))%>%
select(sumry)%>%
unnest()| Gear.Type | Species | effect | group | term | estimate | std.error | statistic |
|---|---|---|---|---|---|---|---|
| GillNet | Hydrophis curtus | fixed | NA | (Intercept) | 0.7987359 | 0.1100043 | 7.260953 |
| GillNet | Hydrophis curtus | fixed | NA | fe | -0.1260171 | 0.0359139 | -3.508871 |
| GillNet | Hydrophis curtus | ran_pars | Boat.Name | sd__(Intercept) | 0.0289787 | NA | NA |
| GillNet | Hydrophis curtus | ran_pars | Residual | sd__Observation | 0.1792341 | NA | NA |
| GillNet | Hydrophis schistosus | fixed | NA | (Intercept) | 0.6309144 | 0.1228337 | 5.136331 |
| GillNet | Hydrophis schistosus | fixed | NA | fe | -0.0596150 | 0.0433414 | -1.375476 |
| GillNet | Hydrophis schistosus | ran_pars | Boat.Name | sd__(Intercept) | 0.0000000 | NA | NA |
| GillNet | Hydrophis schistosus | ran_pars | Residual | sd__Observation | 0.3172430 | NA | NA |
| Trawler | Hydrophis curtus | fixed | NA | (Intercept) | 3.4636588 | 1.7235253 | 2.009636 |
| Trawler | Hydrophis curtus | fixed | NA | fe | -0.3850708 | 0.2993008 | -1.286568 |
| Trawler | Hydrophis curtus | ran_pars | Boat.Name | sd__(Intercept) | 1.3313764 | NA | NA |
| Trawler | Hydrophis curtus | ran_pars | Residual | sd__Observation | 1.1400297 | NA | NA |
| Trawler | Hydrophis schistosus | fixed | NA | (Intercept) | 3.9363791 | 0.9667292 | 4.071853 |
| Trawler | Hydrophis schistosus | fixed | NA | fe | -0.4062698 | 0.1380508 | -2.942899 |
| Trawler | Hydrophis schistosus | ran_pars | Boat.Name | sd__(Intercept) | 0.7191601 | NA | NA |
| Trawler | Hydrophis schistosus | ran_pars | Residual | sd__Observation | 1.4341989 | NA | NA |
CPUE%>%
group_by(Gear.Type, Species)%>%
nest()%>%
mutate(mod = map(data, ~lmer(CPUE ~ fe + (1|Boat.Name), data = .)),
anva = map(mod, car::Anova),
sumry = map(anva, broom::tidy))%>%
select(sumry)%>%
unnest()| Gear.Type | Species | term | statistic | df | p.value |
|---|---|---|---|---|---|
| GillNet | Hydrophis curtus | fe | 12.312177 | 1 | 0.0004500 |
| GillNet | Hydrophis schistosus | fe | 1.891933 | 1 | 0.1689840 |
| Trawler | Hydrophis curtus | fe | 1.655257 | 1 | 0.1982450 |
| Trawler | Hydrophis schistosus | fe | 8.660658 | 1 | 0.0032515 |
Mortality rates of sea snakes in bycatch by gear
We compared risk of mortalities of each species across and within gears using binomial regression with conditions at encounter as the response variable.
#Extracting mortality data from bycatch database
mortality <- snakes%>%
filter(Gear.Type == "GillNet" | Gear.Type == "Trawler")%>%
filter(Species == "Hydrophis curtus" |
Species == "Hydrophis schistosus")%>% # keeping dominant species
group_by(Gear.Type, Species)%>%
summarise(n.alive = sum(Condition.at.encounter. == "Alive"),
n.dead = sum(Condition.at.encounter. == "Dead"),
mort.rate = (n.dead/(n.dead+n.alive))*100)
#table
mortality| Gear.Type | Species | n.alive | n.dead | mort.rate |
|---|---|---|---|---|
| GillNet | Hydrophis curtus | 97 | 34 | 25.954199 |
| GillNet | Hydrophis schistosus | 448 | 34 | 7.053942 |
| Trawler | Hydrophis curtus | 24 | 51 | 68.000000 |
| Trawler | Hydrophis schistosus | 236 | 110 | 31.791908 |
#Plot
mortality%>%
ggplot(aes(Species, mort.rate, fill = Gear.Type))+
geom_col(position = "dodge", col = "black", width = 0.25)+
scale_fill_brewer(palette = "Greys")+
labs(x = "Species", y = "Mortality Rate (%)")+
theme(axis.text.x = element_text(face = "italic"))
ggsave("./Figures/figure 3.tiff", device = "tiff", last_plot(), dpi = "print",
height = 4.5, width = 4.5)# monthly mortalities
mortality.month <- snakes%>%
filter(Gear.Type == "GillNet" | Gear.Type == "Trawler")%>%
filter(Species == "Hydrophis curtus" |
Species == "Hydrophis schistosus")%>% # keeping dominant species
group_by(Month,Gear.Type, Species)%>%
summarise(n.alive = sum(Condition.at.encounter. == "Alive"),
n.dead = sum(Condition.at.encounter. == "Dead"),
mort.rate = (n.dead/(n.dead+n.alive))*100)%>%
ungroup()%>%
mutate(Month = factor(Month, levels = month.abb))%>%
complete(nesting(Gear.Type), Month,
Species, fill = list(mort.rate = 0)) # completing species x trip combinations
mortality.month%>%
mutate(Month = factor(Month, levels = month.abb))%>%
ggplot(aes(Month, mort.rate, shape = Species, group = Species))+
geom_point(size = 3)+
geom_line(aes(linetype = Species), size = 1)+
geom_vline(xintercept = "Jun", size = 1, linetype = "dotted")+
geom_vline(xintercept = "Aug", size = 1, linetype = "dotted")+
geom_text(aes(x = "Jul", y = 50, label = "Monsoon Ban"))+
labs(y = "Mortality Rate (%)")+
facet_wrap(~Gear.Type, ncol = 1)+
theme(legend.text = element_text(face = "italic"))
ggsave("./Figures/figure 3b.tiff", device = "tiff", last_plot(), dpi = "print",
height = 4.5, width = 8)#extracting mortality per trip from bycatch database
mort.grsp <- snakes%>%
filter(Gear.Type == "GillNet" | Gear.Type == "Trawler")%>%
filter(Species == "Hydrophis curtus" |
Species == "Hydrophis schistosus")%>% # keeping dominant species
droplevels()%>%
select(Gear.Type, Date, Boat.Name, Species, Condition.at.encounter.)%>%
rename(Condition = Condition.at.encounter.)%>%
mutate(Condition = as.factor(Condition))# testing difference in mortality between trawlers and gillnets
mort.grsp%>%
group_by(Species)%>%
nest()%>%
mutate(mod = map(data, ~glm(Condition ~ Gear.Type, data = ., family = binomial())),
sumry = map(mod, broom::tidy),
stat = map(mod, ~mcfadden(.)))%>%
select(sumry, stat)%>%
unnest()| Species | term | estimate | std.error | statistic | p.value | stat |
|---|---|---|---|---|---|---|
| Hydrophis schistosus | (Intercept) | -2.578433 | 0.1778870 | -14.494775 | 0e+00 | 0.1131580 |
| Hydrophis schistosus | Gear.TypeTrawler | 1.815081 | 0.2120660 | 8.559038 | 0e+00 | 0.1131580 |
| Hydrophis curtus | (Intercept) | -1.048350 | 0.1993012 | -5.260131 | 1e-07 | 0.1260767 |
| Hydrophis curtus | Gear.TypeTrawler | 1.802122 | 0.3177978 | 5.670657 | 0e+00 | 0.1260767 |
# testing difference in CPUE between Species
mort.grsp%>%
group_by(Gear.Type)%>%
nest()%>%
mutate(mod = map(data, ~glm(Condition ~ Species, data = ., family = binomial())),
sumry = map(mod, broom::tidy),
stat = map(mod, ~mcfadden(.)))%>%
select(sumry, stat)%>%
unnest()| Gear.Type | term | estimate | std.error | statistic | p.value | stat |
|---|---|---|---|---|---|---|
| GillNet | (Intercept) | -1.0483505 | 0.1993014 | -5.260126 | 0.0000001 | 0.0733598 |
| GillNet | SpeciesHydrophis schistosus | -1.5300823 | 0.2671420 | -5.727599 | 0.0000000 | 0.0733598 |
| Trawler | (Intercept) | 0.7537718 | 0.2475368 | 3.045090 | 0.0023261 | 0.0596262 |
| Trawler | SpeciesHydrophis schistosus | -1.5171232 | 0.2731349 | -5.554484 | 0.0000000 | 0.0596262 |
Change in sea snake assemblages
Data on abundance of H. curtus and H. schistosus, overall sampling effort, number of trips was gathered from secondary published; Lobo (2004) and Padate et al (2009). A summary of the number of trips sampled and the sampling effort for each study in Table 2 of supplementary materials A. We used data from trawler bycatch alone to maintain consistency in our comparisons.
#importing data for past studies
past_studies <- read.csv("./Data/Rao et al_data_Past Studies.csv")Sampling effort accross years
#summarising data from current study
effort.studies <- snakes%>%
filter(Gear.Type == "Trawler")%>%
group_by(Year)%>%
distinct(Date, Boat.Name, .keep_all = T)%>%
summarise(N.trips = n(),
Total.fe = sum(Tow.duration.hours, na.rm = T))
#attching to previous studies
effort.studies <- past_studies%>%
select(Year, N.trips, Total.fe)%>%
distinct(Year, .keep_all = T)%>%
bind_rows(effort.studies)
#table
effort.studies| Year | N.trips | Total.fe |
|---|---|---|
| 2002 | 194 | 385.00 |
| 2007 | 70 | 110.00 |
| 2016 | 5 | 34.00 |
| 2017 | 34 | 256.50 |
| 2018 | 26 | 158.66 |
Assemblage composition accross studies
We compared the relative abundance of H. curtus and mortality rates (proportion of individuals encountered dead) across studies using Z-Tests of proportion.
CPUE.studies <- snakes%>%
filter(Gear.Type == "Trawler")%>%
filter(Species == "Hydrophis curtus" |
Species == "Hydrophis schistosus")%>% # keeping dominant species
droplevels()%>%
group_by(Year, Date, Boat.Name, Species)%>%
summarise(Abn = n(), # abundance of sea snakes per trip
fe = mean(Tow.duration.hours, na.rm = T))%>% # fishing effort
complete(nesting(Year, Date, Boat.Name, fe),
Species, fill = list(Abn = 0))%>% # completing species x trip combinations
group_by(Year, Species)%>%
summarise(Mean.CPUE = mean(Abn/fe, na.rm = T), # Mean CPUE per trip
sd = sd(Abn/fe, na.rm = T),
Total.Abn = sum(Abn, na.rm = T),
Total.fe = sum(fe, na.rm = T))
CPUE.studies <- past_studies%>%
select(Year, Species, Mean.CPUE, sd, Total.Abn, Total.fe)%>%
bind_rows(CPUE.studies)#Plotting
assemblage <- CPUE.studies%>%
group_by(Year)%>%
mutate(Total = sum(Total.Abn))%>%
group_by(Year, Species)%>%
summarise(rel.prop = Total.Abn/Total)%>%
ggplot(aes(Year, rel.prop))+
geom_point(aes(shape = Species), size = 3)+
geom_line(aes(linetype = Species, group = Species), size = 1)+
labs(y = "Relative Aundance in byactch")+
theme(axis.text.x = element_blank())+
theme(legend.text = element_text(face = "italic"),
axis.title.x = element_blank())
# adding study periods to plot
study_periods <- read.csv("./Data/Rao et al_data_study_periods.csv")
# rectangle plot
period1 <- ggplot(study_periods)+
geom_rect(aes(xmin = start.year, xmax = end.year, ymax = 1, ymin = -1,
fill = study), col = "black",
size = 0.1, height = 0.1)+
geom_text(aes(x = median.year, y = 0, label = study), size = 3)+
scale_fill_brewer(palette = "Greys", guide = F)+
#scale_y_continuous(limits = c(0.5, 1.3))+
labs(x = "Year")+
theme_void()+
theme(text = element_text(size = 15),
axis.title.x = element_text(size = 15))
# arrow plot
period <- ggplot(study_periods)+
geom_errorbarh(aes(xmin = start.year, xmax = end.year, y = 1.75),
size = 1)+
geom_text(aes(x = median.year, y = 0.75, label = study), size = 3.5)+
scale_y_continuous(limits = c(0, 2.25))+
labs(x = "Year")+
theme_void()+
theme(text = element_text(size =15),
axis.title.x = element_text(size = 15),
axis.text.x = element_text())
# final plot
library(ggpubr)
ggpubr::ggarrange(assemblage, period,
ncol = 1,
heights = c(1, 0.2),
align = "v")
ggsave("./Figures/figure 4.tiff", device = "tiff", last_plot(), dpi = "print",
height = 4.5, width = 8)# bar chart
CPUE.studies%>%
mutate(Year = as.factor(Year))%>%
group_by(Year)%>%
mutate(Total = sum(Total.Abn))%>%
group_by(Year, Species)%>%
summarise(rel.prop = Total.Abn/Total)%>%
ggplot(aes(Year, rel.prop))+
geom_col(aes(fill = Species), col = "black", width = 0.5)+
scale_fill_brewer(palette = "Greys")+
labs(y = "Relative Aundance in byactch")+
theme(legend.text = element_text(face = "italic"))
ggsave("./Figures/figure 4b.tiff", device = "tiff", last_plot(), dpi = "print",
height = 4.5, width = 8)#Z - test for proportion
CPUE.studies%>%
group_by(Year)%>%
mutate(N = sum(Total.Abn))%>%
filter(Species == "Hydrophis curtus")%>%
ungroup()%>%
nest()%>%
mutate(proptest = map(data, ~prop.test(.$Total.Abn, .$N)),
sumry = map(proptest, broom::tidy),
r = map2(proptest, data, ~proptest.r(.x, .y)))%>%
select(sumry,r)%>%
unnest()%>%
select(estimate1:p.value, r)| estimate1 | estimate2 | estimate3 | estimate4 | estimate5 | statistic | p.value | r |
|---|---|---|---|---|---|---|---|
| 0.8495575 | 0.2586207 | 0.1470588 | 0.1764706 | 0.1843318 | 281.6543 | 0 | 0.3995097 |
Difference in mortality of sea snakes in Trawlers across studies
#summarising data from current study
mortality.studies <- snakes%>%
filter(Gear.Type == "Trawler")%>%
filter(Species == "Hydrophis curtus" |
Species == "Hydrophis schistosus")%>% # keeping dominant species
group_by(Year, Species)%>%
summarise(n.alive = sum(Condition.at.encounter. == "Alive"),
n.dead = sum(Condition.at.encounter. == "Dead"),
mort.rate = (n.dead/(n.dead+n.alive)),
Study = "Current")
# Attaching to previous study data
mortality.studies <- past_studies%>%
select(Study, Year, Species, n.alive, n.dead, mort.rate)%>%
bind_rows(mortality.studies)%>%
mutate(N = n.alive + n.dead)%>%
mutate(Study = factor(Study, levels = c("Lobo (2004)", "Padate et al. (2009)", "Current")))%>%
group_by(Study, Species)%>%
summarise(n.dead = sum(n.dead),
N = sum(N),
mort.rate = n.dead/N)
#table
mortality.studies| Study | Species | n.dead | N | mort.rate |
|---|---|---|---|---|
| Lobo (2004) | Hydrophis curtus | 132 | 192 | 0.6875000 |
| Lobo (2004) | Hydrophis schistosus | 4 | 34 | 0.1176471 |
| Padate et al. (2009) | Hydrophis curtus | 15 | NA | NA |
| Padate et al. (2009) | Hydrophis schistosus | 43 | NA | NA |
| Current | Hydrophis curtus | 51 | 75 | 0.6800000 |
| Current | Hydrophis schistosus | 110 | 346 | 0.3179191 |
#plot
mortality.studies%>%
ggplot(aes(Study, mort.rate, fill = Species))+
geom_col(position = position_dodge(preserve = "single"), col = "black", width = 0.25)+
scale_fill_brewer(palette = "Greys")+
labs(x = "Year", y = "Mortality Rate (%)")+
theme(legend.text = element_text(face = "italic"))
# test
mortality.studies%>%
group_by(Species)%>%
drop_na()%>%
nest()%>%
mutate(proptest = map(data, ~prop.test(.$n.dead, .$N)),
#mod = map(data, ~lm(rel.prop ~ Year, data = .)),
sumry = map(proptest, broom::tidy),
r = map2(proptest, data, ~proptest.r(.x, .y)))%>%
select(sumry,r)%>%
unnest()%>%
select(estimate1:p.value, r)| Species | estimate1 | estimate2 | statistic | p.value | r |
|---|---|---|---|---|---|
| Hydrophis curtus | 0.6875000 | 0.6800000 | 0.000000 | 1.0000000 | 0.0000000 |
| Hydrophis schistosus | 0.1176471 | 0.3179191 | 4.997571 | 0.0253829 | 0.0131515 |
Increase in fishing pressure in Konkan 2000 - 2020
In order to understand the change in fishing pressure in the region we gathered data on the number of vessels of each craft type (mechanised, motorised and non- motorised) for Maharashtra and Goa from CMFRI Marine Fisheries Census for 2005 and 2010; Fisheries Department for 2014, 2015, 2017 and 2018; Jena and George (2018) for 2016 and the Handbook on Fisheries Statistics for 2019. The CMFRI Marine Fisheries Census recorded all vessels across multiple harbours across each state. The other sources reported on registered vessels only. Change in the number of vessels of each craft type from 2005 to 2019 were tested using log – linear regression. We compared the number of vessels of each craft across Maharashtra and goa using a t – test.
#importing data on number of vessels
vessels <- read.csv("./Data/Rao et al_data_state-craft-gear.csv")#plot
vessels%>%
ggplot(aes(Year, n, shape = Craft))+
geom_point(size = 3)+
geom_line(aes(group = Craft, linetype = Craft), size = 1)+
scale_y_sqrt()+
labs(y = "No. of Vessels (sq. rt.)")+
facet_wrap(~ State, ncol = 1)
ggsave("./Figures/figure 5.tiff", device = "tiff", last_plot(), dpi = "print",
height = 6, width = 8)#testing change in number of vessels over years
vessels%>%
group_by(State, Craft)%>%
nest()%>%
mutate(mod = map(data, ~glm(n ~ Year, data = ., family = poisson())),
sumry = map(mod, broom::tidy),
r2 = map(mod, ~mcfadden(.)))%>%
select(sumry, r2)%>%
unnest()%>%
arrange(desc(r2))| State | Craft | term | estimate | std.error | statistic | p.value | r2 |
|---|---|---|---|---|---|---|---|
| Maharashtra | Mechanised | (Intercept) | -69.9447237 | 1.3999812 | -49.961190 | 0 | 0.7467118 |
| Maharashtra | Mechanised | Year | 0.0395470 | 0.0006949 | 56.913780 | 0 | 0.7467118 |
| Goa | Mechanised | (Intercept) | -115.4983459 | 4.2006942 | -27.495062 | 0 | 0.5441047 |
| Goa | Mechanised | Year | 0.0611171 | 0.0020846 | 29.318823 | 0 | 0.5441047 |
| Goa | Non - motorised | (Intercept) | 86.2325473 | 8.6498363 | 9.969269 | 0 | 0.3073411 |
| Goa | Non - motorised | Year | -0.0399799 | 0.0042961 | -9.306161 | 0 | 0.3073411 |
| Goa | Motorised | (Intercept) | 252.8613557 | 5.8641018 | 43.120219 | 0 | 0.2632877 |
| Goa | Motorised | Year | -0.1224781 | 0.0029154 | -42.010095 | 0 | 0.2632877 |
| Maharashtra | Motorised | (Intercept) | 185.9700181 | 3.4680508 | 53.623788 | 0 | 0.1044864 |
| Maharashtra | Motorised | Year | -0.0886872 | 0.0017235 | -51.458983 | 0 | 0.1044864 |
| Maharashtra | Non - motorised | (Intercept) | -33.6241328 | 2.2838299 | -14.722696 | 0 | 0.0459403 |
| Maharashtra | Non - motorised | Year | 0.0209820 | 0.0011336 | 18.508824 | 0 | 0.0459403 |
#testing difference in number of vessels between Mah and Goa
vessels%>%
group_by(Craft)%>%
nest()%>%
mutate(mod = map(data, ~t.test(n ~ State, data = .)),
sumry = map(mod, broom::tidy),
d = map(data, ~lsr::cohensD(n ~ State, data = .)))%>%
select(sumry, d)%>%
unnest()%>%
arrange(desc(d))%>%
select(estimate:p.value, d)| Craft | estimate | estimate1 | estimate2 | statistic | p.value | d |
|---|---|---|---|---|---|---|
| Mechanised | -14568.29 | 2054.429 | 16622.714 | -10.598781 | 0.0000233 | 5.6652869 |
| Non - motorised | -5365.00 | 304.625 | 5669.625 | -6.387398 | 0.0003649 | 3.1936991 |
| Motorised | -1103.25 | 558.750 | 1662.000 | -1.260113 | 0.2435439 | 0.6300565 |
Study site map
library(tmap)
library(osmdata)
study_sites <- read.csv("./Data/Rao et al_data_study_site.csv")
# geocoding study sites
library(ggmap)
study_sites <- study_sites%>%
filter(Study != "2006-2008")%>%
mutate(Name = as.character(Name))%>%
mutate_geocode(Name)
study_sites| Study | Name | lon | lat |
|---|---|---|---|
| 2016 - 2018 | Malvan | 73.47105 | 16.06307 |
| 2002 - 2003 | Chapora | 73.74337 | 15.60212 |
| 2002 - 2003 | Vasco da Gama, Mormugao | 73.81133 | 15.39819 |
| 2002 - 2003 | Malim Jetty | 73.83445 | 15.50495 |
| 2002 - 2003 | Betul, Ambelim | 73.95637 | 15.14703 |
| 2016 - 2018 | Vengurla | 73.63890 | 15.85141 |
# converting to simplfe feature
library(sf)
study_sites <- st_as_sf(study_sites, coords = c("lon", "lat"))
st_crs(study_sites) <- 4326
# getting study extent
study_ext <- st_bbox(study_sites)
# getting osm map
library(osmdata)
land_raw <- opq(bbox = study_ext, timeout = 100)%>%
add_osm_feature(key = 'admin_level', value = '4')%>%
osmdata_sf()%>%
unique_osmdata()
land <- land_raw$osm_multipolygons
st_crs(land) <- 4326
india <- read_sf("./Data/India-States.shp")
# plotting site map
ss_main <- tm_shape(land)+
tm_polygons()+
tm_text("name", remove.overlap = F)+
tm_shape(india)+
tm_polygons()+
tm_shape(study_sites, bbox = study_ext+c(-0.5,-.5,0.5,0.5), is.master = T)+
tm_symbols(size = 1, shape = "Study", border.lwd = 2.5, border.col = "black", col = "gray")+
tm_text("Name", size = 0.7,fontface = "bold",
ymod = -0.7, xmod = -0.75, remove.overlap = T)+
tm_graticules(n.x = 5)+
tm_scale_bar(position = c("left", "top"), text.size = 0.5) +
tm_compass(position = c('right', 'top'), size = 5)+
tm_style(style = "gray")+
tm_layout(scale = 1.5, legend.position = c("left", "bottom"))
# making inset map
ext <- st_as_sfc(study_ext+c(-1.2,-1.2,1.2,1.2))
ss_inset <- tm_shape(india)+
tm_polygons()+
tm_shape(ext)+
tm_borders(lwd = 3)+
tm_style("gray")+
tm_layout(scale = 1.5)
# making final map
vp <- grid::viewport(x = 0.75, y = 0.15, width = 0.2, height = 0.3)
tmap_save(ss_main, insets_tm = ss_inset, insets_vp = vp,
"./Figures/Figure_1.tiff", height = 8, width = 8)